library(tidyverse)
library(psych)
library(plotly)
library(GGally)
library(car)
library(gridExtra)
library(cowplot)
library(DescTools)Para a aplicação a seguir, consideraremos o modelo de regressão normal linear
\[ y_i = \beta_1 + \beta_2x_{2i} + ... + \beta_px_{pi} + \epsilon_i, \quad i = 1,...,n \]
em que os erros \(e_i\) são variáveis aleatórias independentes, normalmente distribuídas, de média zero e variância \(\sigma^2\) constante.
Exemplo disponível em: Paula, G. A. (2013). Modelos de Regressão com Apoio Computacional. São Paulo, SP: IME-USP. (Exercício 23, página 112). Dados disponíveis em: <https://www.ime.usp.br/~giapaula/textoregressao.htm>, arquivo imoveis.dat.
Neste exemplo, vamos modelar o preço de venda de imóveis a partir de dados relativos a uma amostra de 27 imóveis. As variáveis do conjunto de dados são:
imposto: imposto do imóvel (em US$ 100);
areaT: área do terreno (em 1.000 pés quadrados);
areaC: área construída (em 1.000 pés quadrados);
idade: idade da residência (em anos);
preco: preço de venda do imóvel (em US$ 1.000).
Sendo assim, o nosso objetivo é encontrar o melhor ajuste linear que explique e quantifique a variável preço de venda a partir das demais.
# Obtendo os dados
dados = read.table("dados/imoveis.dat")
colnames(dados) = c("imposto", "areaT", "areaC", "idade", "preco")
attach(dados)
# Algumas observações dos dados
head(dados)## imposto areaT areaC idade preco
## 1 4.9176 3.472 0.998 42 25.9
## 2 5.0208 3.531 1.500 62 29.5
## 3 4.5429 2.275 1.175 40 27.9
## 4 4.5573 4.050 1.232 54 25.9
## 5 5.0597 4.455 1.121 42 29.9
## 6 3.8910 4.455 0.988 56 29.9
# Medidas descritivas
describe(dados)## vars n mean sd median trimmed mad min max range skew
## imposto 1 27 7.24 2.88 6.09 6.84 2.15 3.89 15.42 11.53 1.40
## areaT 2 27 6.35 2.40 5.85 6.22 2.07 2.28 12.80 10.53 0.68
## areaC 3 27 1.51 0.56 1.49 1.41 0.39 0.98 3.42 2.44 2.02
## idade 4 27 36.48 14.05 40.00 36.96 14.83 3.00 62.00 59.00 -0.34
## preco 5 27 38.50 14.31 36.90 35.65 8.90 25.90 84.90 59.00 2.25
## kurtosis se
## imposto 1.34 0.55
## areaT 0.01 0.46
## areaC 4.08 0.11
## idade -0.54 2.70
## preco 4.63 2.75
plot_ly(type = "box") %>%
add_trace(y = imposto, name = "imposto") %>%
add_trace(y = areaT, name = "área do terreno") %>%
add_trace(y = areaC, name = "área construída") %>%
add_trace(y = idade, name = "idade do imóvel") %>%
add_trace(y = preco, name = "preço de venda")Em um breve resumo descritivo dos dados, observamos que as altas variâncias das variáveis idade e preço destacam-se das demais. Com uma idade mediana de 40 anos, temos um imóvel de apenas 3 anos de idade e, com um preço mediano de US$ 36.900, temos dois imóveis bem mais caros nos valores de US$ 82.900 e US$ 84.900.
Os gráficos de box-plot nos mostram que os dados possuem outliers e que as variáveis possuem distribuições assimétricas. Focando na variável alvo, preço de venda, os pontos destoantes são justamente os dois imóveis mais caros.
g = ggpairs(dados, aes(color = I("slategray"), fill = I("slategray")),
lower = list(continuous = wrap("smooth",col="black")),
diag=list(continuous = wrap("densityDiag",alpha=0.5,size=1)))
ggplotly(g)No gráfico acima, as curvas representadas na diagonal principal deixam mais evidente a constatação de assimetria nas distribuições observada nos box-plots. A variável preço possui correlações positivas fortes com as variáveis imposto, área do terreno e área construída, mas uma correlação negativa fraca com a variável idade. Em outras palavras, quanto maiores os valores de imposto, área construída e área do terreno, maior o preço de venda do imóvel.
Podemos observar também que há correlações relevantes entre as variáveis explicativas imposto, área construída e área do terreno. Isso nos dá indícios de multicolinearidade e, assim, uma possível redundância de informação passada por essas variáveis. Então, talvez um modelo completo com todas as variáveis não seja o mais adequado.
Nesse primeiro ajuste, consideraremos o modelo completo com todas as variáveis disponíveis.
\[ preço_i = \beta_0 + \beta_1imposto_i + \beta_2areaC_i + \beta_3areaT_i + \beta_4idade_i \]
model1 = lm(preco~.,dados)
summary(model1)##
## Call:
## lm(formula = preco ~ ., data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.9685 -2.7152 0.2663 2.1374 6.9872
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.43587 4.09233 0.595 0.55776
## imposto 2.07823 0.55296 3.758 0.00109 **
## areaT 0.23238 0.50658 0.459 0.65094
## areaC 13.97381 2.90663 4.808 8.4e-05 ***
## idade -0.04376 0.06628 -0.660 0.51594
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.077 on 22 degrees of freedom
## Multiple R-squared: 0.9313, Adjusted R-squared: 0.9188
## F-statistic: 74.56 on 4 and 22 DF, p-value: 1.809e-12
Verificamos que, para o modelo completo, as variáveis imposto e área construída são as únicas significativas e obtivemos um \(R^2\) ajustado de 0.92, muito alto, indicando um possível bom ajuste aos dados amostrais. Porém, o \(R^2\) por si só não nos garante que o modelo seja o mais adequado.
Aplicaremos o critério de informação de Akaike (AIC) ao modelo completo, para selecionar as variáveis que deverão permanecer.
MASS::stepAIC(model1)## Start: AIC=80.36
## preco ~ imposto + areaT + areaC + idade
##
## Df Sum of Sq RSS AIC
## - areaT 1 3.50 369.14 78.614
## - idade 1 7.25 372.88 78.887
## <none> 365.64 80.357
## - imposto 1 234.76 600.40 91.748
## - areaC 1 384.13 749.77 97.746
##
## Step: AIC=78.61
## preco ~ imposto + areaC + idade
##
## Df Sum of Sq RSS AIC
## - idade 1 11.51 380.64 77.443
## <none> 369.14 78.614
## - imposto 1 246.10 615.24 90.407
## - areaC 1 489.34 858.47 99.402
##
## Step: AIC=77.44
## preco ~ imposto + areaC
##
## Df Sum of Sq RSS AIC
## <none> 380.64 77.443
## - imposto 1 350.01 730.65 93.049
## - areaC 1 483.16 863.80 97.569
##
## Call:
## lm(formula = preco ~ imposto + areaC, data = dados)
##
## Coefficients:
## (Intercept) imposto areaC
## 0.7903 2.2971 13.9333
O AIC nos retornou as variáveis imposto e área construída como aquelas que deverão ser mantidas no modelo, o qual chamaremos de modelo parcimonioso.
model2 = lm(preco~imposto+areaC,dados)
summary(model2)##
## Call:
## lm(formula = preco ~ imposto + areaC, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.0177 -3.3169 0.0958 2.4382 7.0949
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7903 2.2794 0.347 0.732
## imposto 2.2971 0.4890 4.698 8.96e-05 ***
## areaC 13.9333 2.5244 5.519 1.12e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.982 on 24 degrees of freedom
## Multiple R-squared: 0.9285, Adjusted R-squared: 0.9225
## F-statistic: 155.8 on 2 and 24 DF, p-value: 1.791e-14
E, assim, obtemos um modelo em que as variáveis imposto e área construída são altamente significativas.
Porém, na análise descritiva, observamos que as variáveis área construída e área do terreno possuem uma correlação considerável. Então, podemos ainda testar um segundo modelo parcimonioso considerando a variável área do terreno ao invés da área construída.
model3 = lm(preco~imposto+areaT,dados)
summary(model3)##
## Call:
## lm(formula = preco ~ imposto + areaT, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.621 -3.070 0.446 3.167 10.610
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.1614 3.3019 0.957 0.3479
## imposto 3.9126 0.5293 7.391 1.24e-07 ***
## areaT 1.1017 0.6347 1.736 0.0954 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.655 on 24 degrees of freedom
## Multiple R-squared: 0.8558, Adjusted R-squared: 0.8438
## F-statistic: 71.22 on 2 and 24 DF, p-value: 8.081e-11
E, assim, obtemos um modelo em que a variável área do terreno passa de não significativa para pouco significativa ao nível de 10% de significância. Já a variável imposto continua bastante significativa ao nível de 1% de significância.
Agora, faremos a análise de diagnóstico dos dois modelos parcimoniosos propostos, para verificar suas adequações.
par(mfrow=c(1,2))
qqPlot(model2, pch = 16, main = "preco ~ imposto + areaC")## [1] 6 10
qqPlot(model3, pch = 16, main = "preco ~ imposto + areaT")## [1] 9 27
\(H_0\): Os resíduos têm distribuição normal.
# Shapiro-Wilk
shapiro.test(model2$residuals)##
## Shapiro-Wilk normality test
##
## data: model2$residuals
## W = 0.96558, p-value = 0.4903
shapiro.test(model3$residuals)##
## Shapiro-Wilk normality test
##
## data: model3$residuals
## W = 0.9793, p-value = 0.8462
Ao nível de 5% de significância, não rejeitamos a hipótese de que os resíduos tenham distribuição normal. Nos gráficos de envelope para o ajuste normal, os pontos estão dentro das bandas, exceto pela observação 27 no modelo com área do terreno. Porém, a presença de leves ondulações pode ser indício de variância não constante. Vamos verificar se de fato isso ocorre no gráfico de resíduos por valores ajustados.
residualPlot(model2, type = "rstudent", id = TRUE,
ylab = "resíduo studentizado", xlab = "valor ajustado",
main = "preco ~ imposto + areaC", pch = 16, quadratic = FALSE)residualPlot(model3, type = "rstudent", id = TRUE,
ylab = "resíduo studentizado", xlab = "valor ajustado",
main = "preco ~ imposto + areaT", pch = 16, quadratic = FALSE)Embora não haja violação do pressuposto de normalidade dos resíduos, o gráfico de resíduos por valores ajustados nos dá indícios de heterocedasticidade, isto é, variância não constante, como havíamos suposto a partir do gráfico de envelope. Porém, esse comportamento pode estar sendo influenciado pelas observações destoantes vistas na etapa descritiva. Além disso, podemos perceber no gráfico de resíduos por valores ajustados a presença de pontos aberrantes, isto é, muito além do intervalo [-2,2]. No caso do primeiro modelo parcimonioso, o ponto aberrante é a observação 10, enquanto que no segundo modelo parcimonioso são as observações 9 e 27. Sendo assim, a seguir, vamos verificar se de fato esses pontos estão influenciando na qualidade dos ajustes.
Partiremos agora para a investigação de possíveis pontos influentes, que podem estar atrapalhando a qualidade do ajuste.
Medidas de influência:
DFBetas (dfb): estatísticas que indicam o efeito da remoção de cada observação sobre as estimativas dos parâmetros do modelo;
DFFit (dffit) e Cook’s D (cook.d): são estatísticas que indicam o efeito da remoção de cada observação sobre os valores ajustados do modelo;
COVRATIO (cov.r): estatística que indica o efeito da remoção de cada observação sobre a matriz de covariâncias do modelo, em outras palavras, mede a alteração na precisão das estimativas dos parâmetros do modelo;
HAT (hat): diagonal da matriz de projeção (\(H = X(X'X)^{-1}X'\)) da solução dos mínimos quadrados, é a métrica de alavancagem.
inf = influence.measures(model2)
summary(inf)## Potentially influential observations of
## lm(formula = preco ~ imposto + areaC, data = dados) :
##
## dfb.1_ dfb.imps dfb.areC dffit cov.r cook.d hat
## 9 -0.28 0.00 0.19 0.35 2.18_* 0.04 0.49_*
## 10 -1.25_* 0.29 0.59 1.62_* 0.87 0.73 0.32
## 27 -0.13 -2.04_* 1.86_* -2.12_* 2.04_* 1.39_* 0.61_*
inf = influence.measures(model3)
summary(inf)## Potentially influential observations of
## lm(formula = preco ~ imposto + areaT, data = dados) :
##
## dfb.1_ dfb.imps dfb.areT dffit cov.r cook.d hat
## 9 -1.15_* 1.65_* -0.46 2.00_* 0.81 1.07_* 0.37_*
## 10 -1.21_* 0.47 0.68 1.54_* 1.02 0.69 0.35_*
## 27 0.09 -2.83_* 2.36_* -3.05_* 0.33_* 1.85_* 0.35_*
Baseando-se nas medidas de influência acima, que verificam o efeito da remoção de cada uma das observações individualmente, temos que os possíveis pontos de influência são as observações 9 e 27, para ambos os modelos parcimoniosos, e a observação 10, para o segundo. Percebemos também que a observação 27 é a mais crítica, por impactar nas estimativas dos dois parâmetros, nos valores ajustados e na variância do modelo. Também podemos constatar isso graficamente, como se segue a baixo.
# DISTANCIA DE COOK
plot(model2, which=4, lwd=5, main="Distância de Cook x Observação", caption=F,
sub.caption = "preco ~ imposto + areaC")plot(model3, which=4, lwd=5, main="Distância de Cook x Observação", caption=F,
sub.caption = "preco ~ imposto + areaT")# ALAVANCAGEM
cut = 2*3/27 # 2*p/n
plot(hatvalues(model2), pch = 16,
xlab = "observação", ylab = "medida h", main = "Alavancagem",
sub = "preco ~ imposto + areaC")
abline(h = cut, lty = 2, lwd = 2)
text(hatvalues(model2) * as.integer(hatvalues(model2) > cut),
cex=0.8, p=1, offset=0.3)cut = 2*3/27 # 2*p/n
plot(hatvalues(model3), pch = 16,
xlab = "observação", ylab = "medida h", main = "Alavancagem",
sub = "preco ~ imposto + areaT")
abline(h = cut, lty = 2, lwd = 2)
text(hatvalues(model3) * as.integer(hatvalues(model3) > cut),
cex=0.8, p=1, offset=0.3)As observações 9, 10 e 27 são referentes aos dados:
dados[c(9,10,27),c("imposto","areaC","areaT","preco")]## imposto areaC areaT preco
## 9 15.4202 3.42 9.8 84.9
## 10 14.4598 3.00 12.8 82.9
## 27 12.0000 1.20 5.0 41.0
Os imóveis 9 e 10 possuem os valores de imposto mais elevados e as maiores áreas de terreno e de construção. De acordo com a relação linear identificada na etapa descritiva e a significância dessas variáveis para o modelo preditivo, podemos dizer que essas características justificam um alto preço de venda. Porém, os valores atribuídos aos imóveis 9 e 10 são muito mais elevados que os esperados pelo modelo proposto a imóveis com os mesmos atributos. Já a observação 27 causa estranheza por ter um valor de imposto tão alto e próximo aos dos imóveis 9 e 10, mesmo tendo menos da metade de suas áreas de terreno e de construção e custando pouco menos da metade que o imóvel 10.
# Compara os coeficientes e seus respectivos p-valores de um modelo
# após a remoção de uma observação
impacto = function(m1, m2){
imp = (coef(m1) - coef(m2))/coef(m1) * 100
impacto = data.frame(coef_com_i = coef(m1),
coef_sem_i = coef(m2),
impacto = imp,
pval_com_i = summary(m1)$coefficients[,4],
pval_sem_i = summary(m2)$coefficients[,4])
return(impacto)
}Observação i=9:
model2_sem9 = lm(preco~imposto+areaC,dados,subset=-9)
impacto(model2,model2_sem9)## coef_com_i coef_sem_i impacto pval_com_i pval_sem_i
## (Intercept) 0.7903027 1.433603 -81.39918902 7.318305e-01 0.6304908219
## imposto 2.2970580 2.297773 -0.03113799 8.955810e-05 0.0001221436
## areaC 13.9332510 13.454944 3.43284823 1.122859e-05 0.0001144541
model3_sem9 = lm(preco~imposto+areaT,dados,subset=-9)
impacto(model3,model3_sem9)## coef_com_i coef_sem_i impacto pval_com_i pval_sem_i
## (Intercept) 3.161412 6.558508 -107.45502 3.478823e-01 5.383568e-02
## imposto 3.912561 3.131158 19.97165 1.243498e-07 1.077092e-05
## areaT 1.101699 1.360792 -23.51762 9.543188e-02 2.722222e-02
A remoção da observação 9 não teve impactos relevantes nos coeficientes dos dois modelos parcimoniosos e também não causou mudança inferencial, uma vez que as variáveis imposto, área construída e área do terreno permaneceram igualmente significativas nos respectivos modelos;
Observação i=10:
model2_sem10 = lm(preco~imposto+areaC,dados,subset=-10)
impacto(model2,model2_sem10)## coef_com_i coef_sem_i impacto pval_com_i pval_sem_i
## (Intercept) 0.7903027 3.398847 -330.069021 7.318305e-01 1.640664e-01
## imposto 2.2970580 2.167849 5.624971 8.955810e-05 7.657851e-05
## areaC 13.9332510 12.571459 9.773685 1.122859e-05 2.390378e-05
model3_sem10 = lm(preco~imposto+areaT,dados,subset=-10)
impacto(model3,model3_sem10)## coef_com_i coef_sem_i impacto pval_com_i pval_sem_i
## (Intercept) 3.161412 6.8915170 -117.988556 3.478823e-01 6.463495e-02
## imposto 3.912561 3.6816714 5.901246 1.243498e-07 2.123188e-07
## areaT 1.101699 0.6967438 36.757341 9.543188e-02 2.749152e-01
A remoção da observação 10 não causou impactos relevantes nos coeficientes do primeiro modelo parcimonioso e também não ocasionou mudança inferencial, tendo as variáveis imposto e área construída mantido seus níveis de significância. Já no segundo modelo parcimonioso, a remoção da observação 10 ocasionou uma mudança inferencial, uma vez que a variável área de terreno perde sua significância;
Observação i=27:
model2_sem27 = lm(preco~imposto+areaC,dados,subset=-27)
impacto(model2,model2_sem27)## coef_com_i coef_sem_i impacto pval_com_i pval_sem_i
## (Intercept) 0.7903027 1.068075 -35.14754 7.318305e-01 0.6320276693
## imposto 2.2970580 3.255657 -41.73159 8.955810e-05 0.0001915344
## areaC 13.9332510 9.412139 32.44836 1.122859e-05 0.0155595300
model3_sem27 = lm(preco~imposto+areaT,dados,subset=-27)
impacto(model3,model3_sem27)## coef_com_i coef_sem_i impacto pval_com_i pval_sem_i
## (Intercept) 3.161412 2.9423870 6.928084 3.478823e-01 2.603227e-01
## imposto 3.912561 5.0694489 -29.568550 1.243498e-07 4.765596e-10
## areaT 1.101699 -0.0528502 104.797154 9.543188e-02 9.260592e-01
A remoção da observação 27 impactou em mudança inferencial nos dois modelos parcimoniosos. Enquanto que no primeiro, a significância da variável área construída é reduzida, no segundo, a variável área do terreno perde completamente sua significância.
Como o segundo modelo parcimonioso não nos trouxe nenhuma melhoria para o ajuste e as conclusões sobre os pontos de influência não favoreceram a permanência da variável área do terreno, optaremos pelo primeiro modelo parcimonioso com as variáveis imposto e área construída sugerido inicialmente pelo critério de seleção de Akaike.
A partir dos gráficos de resíduos por valores ajustados, de distância de Cook e de alavancagem, podemos fazer as seguintes classificações:
Ponto aberrante: observação 10;
Pontos de alavanca: observações 9, 10 e 27;
Ponto influente: observação 27.
Então, concluímos que
\[ preço_i = 0.79 + 2.30*imposto_i + 13.93*areaC_i \quad , \]
onde:
A cada US$ 100,00 acrescidos no valor de imposto, aumenta-se, em média, US$ 2.300,00 no preço de venda;
A cada 1.000 metros quadrados de área construída acrescidos, aumenta-se, em média, US$ 13.930,00 no preço de venda.
Exemplo disponível em: Paula, G. A. (2013). Modelos de Regressão com Apoio Computacional. São Paulo, SP: IME-USP. (Exercício 23, página 182). Dados disponíveis em: <https://www.ime.usp.br/~giapaula/textoregressao.htm>, arquivo sobrev.dat.
O conjunto de dados a seguir refere-se à classificação de pacientes com leucemia segundo a ausência ou presença de uma característica morfológica nas células brancas. Os pacientes classificados de AG positivo foram aqueles identificados com a presença dessa característica e os pacientes classificados em AG negativo os que não apresentaram. O conjunto de dados também inclui o tempo de sobrevivência do paciente (em semanas) após o diagnóstico da doença e o número de células brancas (WBC) no momento do diagnóstico.
Variáveis:
WBC: número de células brancas no momento do diagnóstico;
TEMPO: tempo de sobrevivência do paciente (em semanas) após o diagnóstico da doença;
AG: 0 em caso da ausência da característica e 1 em caso da presença.
Supondo que o tempo de sobrevivência após o diagnóstico segue uma distribuição gama, vamos propor um modelo para explicar o tempo médio de sobrevivência dados \(log(WBC)\) e AG (= 1 se positivo, = 0 se negativo).
Denotaremos por \(T_{ij}\) o tempo de sobrevivência do i-ésimo paciente com a j-ésima classificação, \(i = 1,…,33\) e \(j = 1,2\).
# Obtendo os dados
dados = read.csv("dados/sobrev.csv")
dados$AG = factor(dados$AG)
attach(dados)Médias estimadas para o tempo de sobrevivência, considerando o fator AG:
# AG Negativo
dados_ag0 = dados %>% filter(AG==0)
mean(dados_ag0$TEMPO)## [1] 17.9375
# AG Positivo
dados_ag1 = dados %>% filter(AG==1)
mean(dados_ag1$TEMPO)## [1] 62.47059
A partir da diferença entre as médias amostrais acima, podemos supor que os pacientes com o fator AG positivo sobreviveram por mais tempo que aqueles com fator AG negativo.
Coeficientes de variação do tempo de sobrevivência, considerando o fator AG:
# Tempo, AG = 0
sd(dados_ag0$TEMPO)/mean(dados_ag0$TEMPO) * 100## [1] 113.1853
# Tempo, AG = 1
sd(dados_ag1$TEMPO)/mean(dados_ag1$TEMPO) * 100## [1] 87.00598
Ignorando o fator AG, temos o seguinte boxplot:
plot_ly(type = "box", y = TEMPO, name = "TEMPO", color = I("slategray"))Considerando o fator AG, temos os seguintes boxplots:
plot_ly(type = "box",
y = TEMPO, x = AG, name = AG,
color = AG, colors = c("darkred","darkgreen")) %>%
layout(title = "TEMPO x AG")Os gráficos de box-plot evidenciam a assimetria na distribuição da variável tempo de sobrevivência e também apontam a presença de um outlier entre os indivíduos com fator AG negativo. Esse outlier corresponde ao paciente que, mesmo com o fator AG negativo, conseguiu sobreviver a 65 semanas. Nos dados, esse paciente é a 19ª observação.
Ignorando o fator AG, temos a seguinte densidade aproximada do tempo de sobrevivência:
ggplot(dados,aes(x=TEMPO)) +
geom_density(fill = "slategray", color ="slategray",lwd = 1, alpha = 0.5) +
ggtitle("Densidade: TEMPO") +
theme_bw()Considerando o fator AG, temos as seguintes dispersões, densidades e correlações:
g = ggpairs(dados[c("WBC","TEMPO")],
aes(color = AG),
lower = list(continuous = wrap("smooth",col="black")),
diag = list(continuous = wrap("densityDiag", alpha = 0.5, size = 1)))
ggplotly(g)De maneira geral, percebemos uma correlação fraca de -0,33 entre as variáveis WBC e TEMPO. Considerando o fator AG, a correlação é praticamente nula no caso de AG negativo e moderada no caso de AG positivo. Em outras palavras, na presença da característica morfológica, o tempo de sobrevivência tende a diminuir conforme o número de celulas brancas aumentam.
Supondo que \(T_{ij}\) são variáveis aleatórias independentes tais que \(T_{ij} \sim Gama(\mu_{ij},\phi)\) e que \(g(\mu_{ij}) = \eta_{ij}\), com \(\eta_{ij} = \bf{x}_{ij}^\top\beta\), \(\bf{x}_{ij}^\top\) contendo os valores das variáveis explicativas \(log(WBC)\) e \(AG\), e \(\beta\) o vetor de parâmetros.
Assumimos que \(T_{ij}\) segue uma distribuição Gama de média \(\mu_i\) e parâmetro de dispersão \(\phi^{-1}\). Vamos, então, propor um modelo cuja parte sistemática é dada por
\[ \log(\mu_{ij}) = \alpha + \gamma_j + \beta\log(WBC_{ij}) \; , \]
em que \(\gamma_{1} = 0\).
model1 = glm(TEMPO ~ log(WBC) + AG, dados, family = Gamma(link = "log"))
resumo = summary(model1)
resumo##
## Call:
## glm(formula = TEMPO ~ log(WBC) + AG, family = Gamma(link = "log"),
## data = dados)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1746 -1.2663 -0.4251 0.4962 1.9048
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.8154 1.3487 4.312 0.000161 ***
## log(WBC) -0.3044 0.1375 -2.213 0.034617 *
## AG1 1.0176 0.3642 2.794 0.008982 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 1.087715)
##
## Null deviance: 58.138 on 32 degrees of freedom
## Residual deviance: 40.319 on 30 degrees of freedom
## AIC: 301.49
##
## Number of Fisher Scoring iterations: 8
Para o modelo obtido, ao nível de 5% de significância, as variáveis log(WBC) e AG são significativas.
shape = MASS::gamma.shape(model1)
desvio_escal = resumo$deviance * shape$alpha
nivel_descrit = 1 - pchisq(desvio_escal, model1$df.residual)
nivel_descrit## [1] 0.1415408
Também ao nível de 5% de significância, não rejeitamos a hipótese de que o modelo seja adequado.
MASS::stepAIC(model1)## Start: AIC=301.49
## TEMPO ~ log(WBC) + AG
##
## Df Deviance AIC
## <none> 40.319 301.49
## - log(WBC) 1 46.198 304.89
## - AG 1 47.808 306.37
##
## Call: glm(formula = TEMPO ~ log(WBC) + AG, family = Gamma(link = "log"),
## data = dados)
##
## Coefficients:
## (Intercept) log(WBC) AG1
## 5.8154 -0.3044 1.0176
##
## Degrees of Freedom: 32 Total (i.e. Null); 30 Residual
## Null Deviance: 58.14
## Residual Deviance: 40.32 AIC: 301.5
Pelo critério de seleção de Akaike, devemos manter as duas variáveis no modelo.
# ENVELOPE
fit.model = model1
source("envel_gama")Pelo gráfico de envelope, verificamos que todos os pontos se encontram dentro das bandas.
par(mfrow=c(2,2))
# RESÍDUOS STUDENTIZADOS
residualPlot(model1, type = "deviance", id = TRUE,
ylab = "res. compon. do desvio", xlab = "valor ajustado",
pch = 16, smooth = F)
abline(h = c(-2,2), col="red")
# VARIÁVEL Z
w = model1$weights
eta = model1$linear.predictors
z = eta + resid(model1, type="pearson")/sqrt(w)
plot(model1$linear.predictors, z,
xlab = "Preditor Linear", ylab = "Variável z", pch = 16)
# ALAVANCAGEM
cut = min(sort(hatvalues(model1), decreasing = TRUE)[1:3])
plot(hatvalues(model1), pch = 16,
xlab = "observação", ylab = "medida h", main = "Alavancagem")
text(hatvalues(model1) * as.integer(hatvalues(model1) >= cut),
cex=0.8, p=1, offset=0.3)
# DISTÂNCIA DE COOK
plot(model1, which=4, lwd=3, main = "Distância de Cook", caption = F)No gráfico de resíduos por valores ajustados, percebemos uma boa distribuição dos pontos no intervalo de -2 a 2. Embora seja perceptível também um afunilamento à direita, a quantidade de pontos nessa região não é suficiente para indicar que seja uma tendência. Não há presença de pontos aberrantes, uma vez que todos os pontos estão dentro ou muito próximo do intervalo de -2 a 2. A observação 33 fica em evidência no gráfico de distância de Cook, mas não a consideraremos para análise uma vez que a distância não é expressiva em escala. No gráfico de alavancagem, a observação 2 se destaca e, sendo assim, investigaremos sua possível influência na qualidade do ajuste.
summary(influence.measures(model1))## Potentially influential observations of
## glm(formula = TEMPO ~ log(WBC) + AG, family = Gamma(link = "log"), data = dados) :
##
## dfb.1_ dfb.l(WB dfb.AG1 dffit cov.r cook.d hat
## 2 0.09 -0.09 0.04 0.11 1.37_* 0.01 0.20
As medidas de influência acima apontam apenas para a observação 2 como possível ponto influente.
model1_sem2 = glm(TEMPO ~ log(WBC) + AG, dados, subset = -2,
family = Gamma(link = "log"))
impacto(model1, model1_sem2)## coef_com_i coef_sem_i impacto pval_com_i pval_sem_i
## (Intercept) 5.8154109 5.6913329 2.1336073 0.0001610584 0.000584179
## log(WBC) -0.3044004 -0.2919183 4.1005537 0.0346169527 0.062619407
## AG1 1.0176452 1.0090587 0.8437636 0.0089824909 0.011413245
Ao nível de 5% de significância, verificamos que a variável log(WBC) deixa de ser significativa quando removemos a observação 2. Sendo assim, ela causa mudança infefrencial e classifica-se como ponto de alavanca e influente.
Com essa análise, chegamos ao seguinte modelo para explicar o tempo médio de sobrevivência:
\[ \begin{aligned} \log(\mu_{i1}) &= 5.8 - 0.3\log(WBC_{i1}) \; ; \\ \log(\mu_{i2}) &= 6.8 - 0.3\log(WBC_{i2}) \; . \end{aligned} \]
Fixando-se o fator \(\gamma_{j}\), temos que
\[ \frac{\mu_{j}(x+1)}{\mu_{j}(x)} = \frac{\exp{(5.8+ \gamma_j -0.3(x+1))}}{\exp{(5.8 + \gamma_j - 0.3x)}} =\exp{(-0.3) = 0.74} \; . \]
Alterando-se o fator \(\gamma_j\), temos que
\[ \frac{\mu_{i2}}{\mu_{i1}} = \frac{\exp{(5.8 + \gamma_2 - 0.3X_{i2})}}{\exp{(5.8 + \gamma_1 - 0.3X_{i1})}} = \exp{(\gamma_2)} = 2.72 \]
De onde concluímos que, ao fixarmos o fator característica morfológica, o acréscimo de 1 unidade à variável \(\log(WBC)\) representa, em média, uma redução de 26% do tempo de sobrevivência. Já considerando a característica morfológica, o fato do paciente apresentá-la representa um aumento médio de 172% no tempo de sobrevivência.
Os conjuntos de dados apresentados a seguir são provenientes de um experimento dose-resposta, conduzido para avaliar a influência dos extratos vegetais “aquoso frio de folhas”, “aquoso frio de frutos” e um extrato químico, respectivamente, na morte de um determinado tipo de caramujo. Para cada conjunto, ajustaremos um modelo logístico linear simples e um modelo log-log linear simples. Para o melhor ajuste, usando envelope como critério, encontraremos um intervalo assintótico de 95% para a dose letal \(DL_{50}\).
Dados do Extrato Aquoso Frio de Folhas:
dose1 = read.table("dados/dose1.dat", col.names = c("dose","expostos","mortos"))
dose1 = dose1 %>% mutate(restantes = expostos - mortos,
prop.mort = mortos/expostos,
prop.rest = restantes/expostos)
dose1## dose expostos mortos restantes prop.mort prop.rest
## 1 0 50 4 46 0.08 0.92
## 2 15 50 5 45 0.10 0.90
## 3 20 50 14 36 0.28 0.72
## 4 25 50 29 21 0.58 0.42
## 5 30 50 38 12 0.76 0.24
## 6 35 50 41 9 0.82 0.18
## 7 40 50 47 3 0.94 0.06
Dados do Extrato Aquoso Frio de Frutos:
dose2 = read.table("dados/dose2.dat", col.names = c("dose","expostos","mortos"))
dose2 = dose2 %>% mutate(restantes = expostos - mortos,
prop.mort = mortos/expostos,
prop.rest = restantes/expostos)
dose2## dose expostos mortos restantes prop.mort prop.rest
## 1 0 65 2 63 0.03076923 0.96923077
## 2 100 52 2 50 0.03846154 0.96153846
## 3 150 52 8 44 0.15384615 0.84615385
## 4 200 52 20 32 0.38461538 0.61538462
## 5 250 52 41 11 0.78846154 0.21153846
## 6 300 52 48 4 0.92307692 0.07692308
## 7 350 52 51 1 0.98076923 0.01923077
## 8 400 52 52 0 1.00000000 0.00000000
Dados do Extrato Químico:
dose3 = read.table("dados/dose3.dat", col.names = c("dose","expostos","mortos"))
dose3 = dose3 %>% mutate(restantes = expostos - mortos,
prop.mort = mortos/expostos,
prop.rest = restantes/expostos)
dose3## dose expostos mortos restantes prop.mort prop.rest
## 1 0.000 30 2 28 0.06666667 0.9333333
## 2 0.010 30 4 26 0.13333333 0.8666667
## 3 0.015 30 4 26 0.13333333 0.8666667
## 4 0.020 30 11 19 0.36666667 0.6333333
## 5 0.030 29 9 20 0.31034483 0.6896552
## 6 0.040 30 14 16 0.46666667 0.5333333
## 7 0.050 30 24 6 0.80000000 0.2000000
p1 = ggplot(dose1, aes(x = dose, y = prop.mort)) +
geom_col(fill = "darkred") +
ylab("prop. mortos") +
labs(title = "Extrato Aquoso Frio de Folhas") +
theme_bw()
p2 = ggplot(dose2, aes(x = dose, y = prop.mort)) +
geom_col(fill = "darkred") +
ylab("prop. mortos") +
labs(title = "Extrato Aquoso Frio de Frutos") +
theme_bw()
p3 = ggplot(dose3, aes(x = dose, y = prop.mort)) +
geom_col(fill = "darkred") +
ylab("prop. mortos") +
labs(title = "Extrato Químico") +
theme_bw()
plot_grid(p1,p2,p3,ncol=2)# EAF Folhas
log1 = glm(as.matrix(dose1[3:4]) ~ dose, data = dose1, family=binomial)
summary(log1)##
## Call:
## glm(formula = as.matrix(dose1[3:4]) ~ dose, family = binomial,
## data = dose1)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7
## 2.1850 -1.7488 -0.9036 0.7107 0.7624 -0.4650 0.4818
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.80668 0.45190 -8.424 <2e-16 ***
## dose 0.15707 0.01707 9.199 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 172.446 on 6 degrees of freedom
## Residual deviance: 10.184 on 5 degrees of freedom
## AIC: 40.101
##
## Number of Fisher Scoring iterations: 5
# EAF Frutos
log2 = glm(as.matrix(dose2[3:4]) ~ dose, data = dose2, family=binomial)
summary(log2)##
## Call:
## glm(formula = as.matrix(dose2[3:4]) ~ dose, family = binomial,
## data = dose2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.8190 -0.3956 0.1640 0.7142 2.1973
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.546200 0.565304 -9.811 <2e-16 ***
## dose 0.026539 0.002561 10.363 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 353.3692 on 7 degrees of freedom
## Residual deviance: 7.0361 on 6 degrees of freedom
## AIC: 33.504
##
## Number of Fisher Scoring iterations: 5
# E Químico
log3 = glm(as.matrix(dose3[3:4]) ~ dose, data = dose3, family=binomial)
summary(log3)##
## Call:
## glm(formula = as.matrix(dose3[3:4]) ~ dose, family = binomial,
## data = dose3)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7
## -0.05150 0.03496 -0.65291 1.61394 -0.83739 -1.02803 0.98908
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.6016 0.3662 -7.105 1.2e-12 ***
## dose 71.0948 11.3952 6.239 4.4e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 55.0428 on 6 degrees of freedom
## Residual deviance: 5.7713 on 5 degrees of freedom
## AIC: 33.346
##
## Number of Fisher Scoring iterations: 4
# EAF Folhas
1 - pchisq(log1$deviance, log1$df.residual)## [1] 0.07019215
# EAF Frutos
1 - pchisq(log2$deviance, log1$df.residual)## [1] 0.2179659
# E Químico
1 - pchisq(log3$deviance, log1$df.residual)## [1] 0.3291096
# EAF Folhas
cloglog1 = glm(as.matrix(dose1[3:4]) ~ dose, data = dose1,
family = binomial(link = "cloglog"))
summary(cloglog1)##
## Call:
## glm(formula = as.matrix(dose1[3:4]) ~ dose, family = binomial(link = "cloglog"),
## data = dose1)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7
## 0.9330 -2.0977 -0.7403 1.3129 1.3331 -0.5276 -0.5728
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.99021 0.30754 -9.723 <2e-16 ***
## dose 0.10348 0.01017 10.175 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 172.4463 on 6 degrees of freedom
## Residual deviance: 9.9264 on 5 degrees of freedom
## AIC: 39.844
##
## Number of Fisher Scoring iterations: 4
# EAF Frutos
cloglog2 = glm(as.matrix(dose2[3:4]) ~ dose, data = dose2,
family = binomial(link = "cloglog"))
summary(cloglog2)##
## Call:
## glm(formula = as.matrix(dose2[3:4]) ~ dose, family = binomial(link = "cloglog"),
## data = dose2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.48238 -0.94709 -0.03929 0.18022 1.80452
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.988340 0.372324 -10.71 <2e-16 ***
## dose 0.016428 0.001498 10.96 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 353.3692 on 7 degrees of freedom
## Residual deviance: 8.6014 on 6 degrees of freedom
## AIC: 35.069
##
## Number of Fisher Scoring iterations: 7
# E Químico
cloglog3 = glm(as.matrix(dose3[3:4]) ~ dose, data = dose3,
family = binomial(link = "cloglog"))
summary(cloglog3)##
## Call:
## glm(formula = as.matrix(dose3[3:4]) ~ dose, family = binomial(link = "cloglog"),
## data = dose3)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7
## -0.25501 -0.02679 -0.62490 1.74222 -0.60022 -0.90515 0.58896
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.4987 0.3128 -7.988 1.38e-15 ***
## dose 56.7822 8.4721 6.702 2.05e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 55.043 on 6 degrees of freedom
## Residual deviance: 5.018 on 5 degrees of freedom
## AIC: 32.593
##
## Number of Fisher Scoring iterations: 4
# EAF Folhas
1 - pchisq(cloglog1$deviance, cloglog1$df.residual)## [1] 0.07734658
# EAF Frutos
1 - pchisq(cloglog2$deviance, cloglog2$df.residual)## [1] 0.1972691
# E químico
1 - pchisq(cloglog3$deviance, cloglog3$df.residual)## [1] 0.4136878
Extrato Aquoso Frio de Folhas:
fit.model = log1
ntot = dose1$expostos
source("envelr_bino.txt")fit.model = cloglog1
ntot = dose1$expostos
source("envelr_bino_cloglog.txt")Extrato Aquoso Frio de Frutos:
fit.model = log2
ntot = dose2$expostos
source("envelr_bino.txt")fit.model = cloglog2
ntot = dose2$expostos
source("envelr_bino_cloglog.txt")Extrato Químico:
fit.model = log3
ntot = dose3$expostos
source("envelr_bino.txt")fit.model = cloglog3
ntot = dose3$expostos
source("envelr_bino_cloglog.txt")pi_log1 = exp(log1$coeff[1] + log1$coeff[2] * dose1$dose) /
(1 + exp(log1$coeff[1] + log1$coeff[2] * dose1$dose))
p_log1 = ggplot(dose1, aes(x = dose, y = prop.mort)) +
geom_point() +
geom_line(y = pi_log1) +
ylab("prop. mortos") +
ggtitle("E.A.F.Folhas - Log simples") +
theme_bw()
pi_cloglog1 = exp(cloglog1$coeff[1] + cloglog1$coeff[2]*dose1$dose) /
(1 + exp(cloglog1$coeff[1] + cloglog1$coeff[2]*dose1$dose))
p_cloglog1 = ggplot(dose1, aes(x = dose, y = prop.mort)) +
geom_point() +
geom_line(y = pi_cloglog1) +
ylab("prop. mortos") +
ggtitle("E.A.F.Folhas - Complem. Log-log") +
theme_bw()
pi_log2 = exp(log2$coeff[1] + log2$coeff[2] * dose2$dose) /
(1 + exp(log2$coeff[1] + log2$coeff[2] * dose2$dose))
p_log2 = ggplot(dose2, aes(x = dose, y = prop.mort)) +
geom_point() +
geom_line(y = pi_log2) +
ylab("prop. mortos") +
ggtitle("E.A.F.Frutos - Log simples") +
theme_bw()
pi_cloglog2 = exp(cloglog2$coeff[1] + cloglog2$coeff[2]*dose2$dose) /
(1 + exp(cloglog2$coeff[1] + cloglog2$coeff[2]*dose2$dose))
p_cloglog2 = ggplot(dose2, aes(x = dose, y = prop.mort)) +
geom_point() +
geom_line(y = pi_cloglog2) +
ylab("prop. mortos") +
ggtitle("E.A.F.Frutos - Complem. Log-log") +
theme_bw()
pi_log3 = exp(log3$coeff[1] + log3$coeff[2] * dose3$dose) /
(1 + exp(log3$coeff[1] + log3$coeff[2] * dose3$dose))
p_log3 = ggplot(dose3, aes(x = dose, y = prop.mort)) +
geom_point() +
geom_line(y = pi_log3) +
ylab("prop. mortos") +
ggtitle("E. Químico - Log simples") +
theme_bw()
pi_cloglog3 = exp(cloglog3$coeff[1] + cloglog3$coeff[2]*dose3$dose) /
(1 + exp(cloglog3$coeff[1] + cloglog3$coeff[2]*dose3$dose))
p_cloglog3 = ggplot(dose3, aes(x = dose, y = prop.mort)) +
geom_point() +
geom_line(y = pi_cloglog3) +
ylab("prop. mortos") +
ggtitle("E. Químico - Complem. Log-log") +
theme_bw()
plot_grid(p_log1,p_cloglog1,
p_log2,p_cloglog2,
p_log3,p_cloglog3, ncol = 2)# DOSE LETAL ESTIMATIVA PONTUAL
DL_100p = function(model,p){
(1/model$coeff[2])*(log(p/(1-p)) - model$coeff[1])
}
# DOSE LETAL ESTIMATIVA INTERVALAR, 95% DE CONFIANÇA
DL_IC_log = function(model, p, gama){
db = cbind(-1/model$coeff[2],
1/(model$coeff[2]^2)*(model$coeff[1]-log(p/(1-p))))
X = model.matrix(model)
w = model$weights
W = diag(w)
Var = db %*% solve(t(X) %*% W %*% X) %*% t(db)
LI = DL_100p(model, 0.5) - qnorm((1+gama)/2)*sqrt(Var)
LS = DL_100p(model, 0.5) + qnorm((1+gama)/2)*sqrt(Var)
return(c(LI,LS))
}data.frame(extrato = c("EAF de Folhas","EAF de Frutos", "E Químico"),
LI = c(DL_IC_log(log1, 0.5, 0.95)[1],
DL_IC_log(log2, 0.5, 0.95)[1],
DL_IC_log(log3, 0.5, 0.95)[1]),
DL_100p = c(DL_100p(log1,0.5),DL_100p(log2,0.5),DL_100p(log3,0.5)),
LS = c(DL_IC_log(log1, 0.5, 0.95)[2],
DL_IC_log(log2, 0.5, 0.95)[2],
DL_IC_log(log3, 0.5, 0.95)[2]))## extrato LI DL_100p LS
## 1 EAF de Folhas 22.50628769 24.23566646 25.96504522
## 2 EAF de Frutos 197.03961411 208.98193852 220.92426293
## 3 E Químico 0.03131498 0.03659294 0.04187089
Exemplo disponível em: Paula, G. A. (2013). Modelos de Regressão com Apoio Computacional. São Paulo, SP: IME-USP. (Exercício 20, página 277). Dados disponíveis em: <https://www.ime.usp.br/~giapaula/textoregressao.htm>, arquivo leuc.dat.
Os dados a seguir referem-se a um estudo com 51 pacientes adultos, previamente diagnosticados com um tipo agudo de leucemia, e que receberam um tipo de tratamento. A eficiência ou não do tratamento foi verificada após um certo período. As variáveis em estudo são:
idade (\(x_1\)): idade do paciente (em anos);
mancha (\(x_2\)): mancha diferencial da doença (em %);
infilt (\(x_3\)): infiltração na medula (em %);
leuc (\(x_4\)): células com leucemia (em %);
malig (\(x_5\)): malignidade da doença (\(\times 10^3\));
tmax (\(x_6\)): temperatura máxima pré-tratamento (\(\times 10^\circ F\));
trat (\(y\)): tratamento (0: não-satisfatório, 1: satisfatório).
Consideremos o seguinte modelo logístico linear, para explicar a probabilidade de eficiência do tratamento a partir das seis primeiras variáveis explicativas acima:
\[ Y \sim Bernoulli(\pi(\mathbf{x})) \]
\[ \log \left\{ \frac{\pi(\mathbf{x})}{1-\pi(\mathbf{x})} \right\} = \beta_1 + \beta_2x_2 + ... + \beta_6x_6 \]
dados = read.table("dados/leuce.dat",
col.names = c("idade","mancha","infilt",
"leuc","malig","tmax",
"trat","tvida","status"))
dados$trat = factor(dados$trat)
dados = dados[1:7]
attach(dados)p1 = dados %>% ggplot(aes(x = trat, y = idade, fill = trat)) +
geom_boxplot() +
scale_fill_manual(values=c("#BC3C29", "#20854E")) +
ggtitle("idade x tratamento") +
theme(legend.position = "none")
p2 = dados %>% ggplot(aes(x = trat, y = mancha, fill = trat)) +
geom_boxplot() +
scale_fill_manual(values=c("#BC3C29", "#20854E")) +
ggtitle("mancha diferencial x tratamento") +
theme(legend.position = "none")
p3 = dados %>% ggplot(aes(x = trat, y = infilt, fill = trat)) +
geom_boxplot() +
scale_fill_manual(values=c("#BC3C29", "#20854E")) +
ggtitle("infiltração x tratamento") +
theme(legend.position = "none")
p4 = dados %>% ggplot(aes(x = trat, y = leuc, fill = trat)) +
geom_boxplot() +
scale_fill_manual(values=c("#BC3C29", "#20854E")) +
ggtitle("células leucêmicas x tratamento") +
theme(legend.position = "none")
p5 = dados %>% ggplot(aes(x = trat, y = malig, fill = trat)) +
geom_boxplot() +
scale_fill_manual(values=c("#BC3C29", "#20854E")) +
ggtitle("malignidade x tratamento") +
theme(legend.position = "none")
p6 = dados %>% ggplot(aes(x = trat, y = tmax, fill = trat)) +
geom_boxplot() +
scale_fill_manual(values=c("#BC3C29", "#20854E")) +
ggtitle("temperatura máxima x tratamento") +
theme(legend.position = "none")
grid.arrange(p1,p2,p3,p4,p5,p6)Observando os box-plots acima, fica evidente que o tratamento foi mais eficaz em pacientes mais jovens e naqueles com maior porcentagem de células leucêmicas.
corrplot::corrplot(cor(dados[1:6]), method = "color", type = "lower",
addCoef.col = "white", tl.col = "black")Dentre as variáveis explicativas, o par mancha diferencial e infiltração na medula foi o único com alto grau de correlação. Então, vale estudar modelos considerando as duas individualmente junto as demais.
model = glm(trat ~ idade + mancha + infilt + leuc + malig + tmax,
family = binomial)
summary(model)##
## Call:
## glm(formula = trat ~ idade + mancha + infilt + leuc + malig +
## tmax, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.56756 -0.55762 -0.05269 0.59038 2.38751
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 108.33115 41.84379 2.589 0.00963 **
## idade -0.06231 0.02746 -2.269 0.02327 *
## mancha -0.00469 0.04005 -0.117 0.90677
## infilt 0.03104 0.03789 0.819 0.41264
## leuc 0.37281 0.13247 2.814 0.00489 **
## malig 0.03267 0.04605 0.710 0.47801
## tmax -0.11162 0.04263 -2.618 0.00884 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 70.524 on 50 degrees of freedom
## Residual deviance: 39.275 on 44 degrees of freedom
## AIC: 53.275
##
## Number of Fisher Scoring iterations: 6
Para o modelo com todas as variáveis fornecidas, mancha, infiltração e malignidade não foram significativas.
MASS::stepAIC(model)## Start: AIC=53.28
## trat ~ idade + mancha + infilt + leuc + malig + tmax
##
## Df Deviance AIC
## - mancha 1 39.289 51.289
## - infilt 1 40.014 52.014
## - malig 1 40.115 52.115
## <none> 39.275 53.275
## - idade 1 45.737 57.737
## - tmax 1 48.760 60.760
## - leuc 1 52.213 64.213
##
## Step: AIC=51.29
## trat ~ idade + infilt + leuc + malig + tmax
##
## Df Deviance AIC
## - malig 1 40.136 50.136
## - infilt 1 41.042 51.042
## <none> 39.289 51.289
## - idade 1 45.784 55.784
## - tmax 1 48.882 58.882
## - leuc 1 52.728 62.728
##
## Step: AIC=50.14
## trat ~ idade + infilt + leuc + tmax
##
## Df Deviance AIC
## <none> 40.136 50.136
## - infilt 1 43.265 51.265
## - idade 1 46.438 54.438
## - tmax 1 48.971 56.971
## - leuc 1 57.602 65.602
##
## Call: glm(formula = trat ~ idade + infilt + leuc + tmax, family = binomial)
##
## Coefficients:
## (Intercept) idade infilt leuc tmax
## 95.56766 -0.06026 0.03413 0.40673 -0.09944
##
## Degrees of Freedom: 50 Total (i.e. Null); 46 Residual
## Null Deviance: 70.52
## Residual Deviance: 40.14 AIC: 50.14
No método stepwise AIC, as variáveis selecionadas foram idade, infilt, leuc e tmax. Como infilt e mancha são fortemente correlacionadas, vamos testar também um modelo com mancha ao invés de infilt.
# Com a variável INFILT
modelr1 = glm(trat ~ idade + infilt + leuc + tmax, family = binomial)
summary(modelr1)##
## Call:
## glm(formula = trat ~ idade + infilt + leuc + tmax, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.73886 -0.56473 -0.05442 0.62185 2.26516
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 95.56766 38.59482 2.476 0.01328 *
## idade -0.06026 0.02678 -2.250 0.02445 *
## infilt 0.03413 0.02079 1.641 0.10077
## leuc 0.40673 0.13034 3.121 0.00181 **
## tmax -0.09944 0.03954 -2.515 0.01191 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 70.524 on 50 degrees of freedom
## Residual deviance: 40.136 on 46 degrees of freedom
## AIC: 50.136
##
## Number of Fisher Scoring iterations: 6
# Com a variável MANCHA
modelr2 = glm(trat ~ idade + mancha + leuc + tmax, family = binomial)
summary(modelr2)##
## Call:
## glm(formula = trat ~ idade + mancha + leuc + tmax, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.87043 -0.58810 -0.05384 0.65397 2.30025
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 94.45927 37.25941 2.535 0.01124 *
## idade -0.05637 0.02632 -2.141 0.03225 *
## mancha 0.03029 0.02191 1.382 0.16682
## leuc 0.41243 0.12938 3.188 0.00143 **
## tmax -0.09854 0.03826 -2.576 0.01000 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 70.524 on 50 degrees of freedom
## Residual deviance: 41.222 on 46 degrees of freedom
## AIC: 51.222
##
## Number of Fisher Scoring iterations: 6
Os modelo reduzido com a variável infilt quanto o com a variável mancha apresentam as mesmas estimativas aproximadas para os parâmetros. Também nos dois, essas variáveis permanecem não significativas. Adiante, estudaremos uma possível influência de pontos sobre a inferência do modelo e avaliaremos a remoção das variáveis.
1 - pchisq(modelr1$deviance, modelr1$df.residual)## [1] 0.7153545
1 - pchisq(modelr2$deviance, modelr2$df.residual)## [1] 0.6723047
Ao nível de 5% de significância, não rejeitamos a hipótese de que os modelos sejam adequados.
fit.model = modelr1
source("envel_bino.txt")fit.model = modelr2
source("envel_bino.txt")Pelos gráficos de envelope, temos que o modelo reduzido com a variável mancha ficou melhor ajustado e, então, o escolheremos para prosseguir nossa análise.
par(mfrow=c(2,2))
# ALAVANCAGEM
cut = min(sort(hatvalues(modelr2), decreasing = TRUE)[1:3])
plot(hatvalues(modelr2), pch = 16,
xlab = "observação", ylab = "medida h", main = "Alavancagem")
text(hatvalues(modelr2) * as.integer(hatvalues(modelr2) >= cut),
cex=0.8, p=1, offset=0.4)
# DISTÂNCIA DE COOK
plot(modelr2, which=4, lwd=3, main = "Distância de Cook", caption = F)
# RESÍDUO COMPONENTE DO DESVIO POR ÍNDICE
plot(residuals(modelr2), pch = 16,
ylab = "res. compon. do desvio", xlab = "índice")
abline(h = c(-2,2), col="red")No gráfico de resíduo componente do desvio, verificamos que os pontos estão bem distribuídos aleatoriamente e sem nenhuma observação aberrante muito além dos limites de confiança. Nos gráficos de alavancagem e distância de Cook, destacam–se os pontos 2, 13, 32, 47 e 48. A seguir, verificaremos se há influência desses pontos sobre a inferência do modelo.
summary(influence.measures(modelr2))## Potentially influential observations of
## glm(formula = trat ~ idade + mancha + leuc + tmax, family = binomial) :
##
## dfb.1_ dfb.idad dfb.mnch dfb.leuc dfb.tmax dffit cov.r cook.d hat
## 2 -0.61 -0.36 -0.22 0.02 0.62 1.06_* 1.31 0.14 0.33_*
## 13 -0.09 0.23 0.60 0.14 0.05 -0.74 1.33_* 0.06 0.28
## 32 0.23 -0.16 -0.44 -0.46 -0.19 -1.03_* 0.87 0.20 0.20
## 47 -0.45 0.58 -0.22 -0.40 0.44 0.71 0.57_* 0.20 0.07
## 48 -0.30 -0.34 0.55 -0.44 0.29 -1.11_* 1.11 0.18 0.28
Observação 2
modelr2_sem2 = glm(trat ~ idade + mancha + leuc + tmax, subset = -2,
family = binomial)
impacto(modelr2, modelr2_sem2)[3:5]## impacto pval_com_i pval_sem_i
## (Intercept) -19.7958316 0.011238999 0.007759372
## idade 12.2404843 0.032251493 0.061370014
## mancha -13.4866607 0.166824037 0.127772519
## leuc -0.7081507 0.001434345 0.001634411
## tmax -19.8367890 0.010004635 0.007094660
A 5% de significância, removendo-se a observação 2, a variável idade deixa de ser significativa.
Observação 13
modelr2_sem13 = glm(trat ~ idade + mancha + leuc + tmax, subset = -13,
family = binomial)
impacto(modelr2, modelr2_sem13)[3:5]## impacto pval_com_i pval_sem_i
## (Intercept) -3.204997 0.011238999 0.009617055
## idade -8.835451 0.032251493 0.024988565
## mancha 32.179710 0.166824037 0.393890235
## leuc 2.477318 0.001434345 0.001728282
## tmax -1.989041 0.010004635 0.008939169
A 5% de significância, removendo-se a observação 13, não há mudança inferencial no modelo.
Observação 32
modelr2_sem32 = glm(trat ~ idade + mancha + leuc + tmax, subset = -32,
family = binomial)
impacto(modelr2, modelr2_sem32)[3:5]## impacto pval_com_i pval_sem_i
## (Intercept) 8.654518 0.011238999 0.026129747
## idade 7.325326 0.032251493 0.052668960
## mancha -36.997935 0.166824037 0.091633885
## leuc -17.602426 0.001434345 0.001022189
## tmax 6.768428 0.010004635 0.021156787
A 5% de significância, removendo-se a observação 32, a variável idade deixa de ser significativa.
Observação 47
modelr2_sem47 = glm(trat ~ idade + mancha + leuc + tmax, subset = -47,
family = binomial)
impacto(modelr2, modelr2_sem47)[3:5]## impacto pval_com_i pval_sem_i
## (Intercept) -33.67437 0.011238999 0.006535694
## idade -50.54577 0.032251493 0.013157635
## mancha -29.65272 0.166824037 0.103656074
## leuc -26.05928 0.001434345 0.001848562
## tmax -33.01815 0.010004635 0.006022423
A 5% de significância, removendo-se a observação 47, não há mudança inferencial no modelo.
Observação 48
modelr2_sem48 = glm(trat ~ idade + mancha + leuc + tmax, subset = -48,
family = binomial)
impacto(modelr2, modelr2_sem48)[3:5]## impacto pval_com_i pval_sem_i
## (Intercept) -13.22389 0.011238999 0.009203440
## idade 12.27356 0.032251493 0.061894967
## mancha 33.26921 0.166824037 0.371617432
## leuc -14.88983 0.001434345 0.001479636
## tmax -12.88917 0.010004635 0.008270069
A 5% de significância, removendo-se a observação 48, a variável idade deixa de ser significativa.
As variáveis mancha e infilt não foram significativas em nenhuma circustância vista anteriormente. Então, podemos considerar removê-las. Além disso, a variável idade deixou de ser significativa com a remoção de algumas observações. Então, testaremos a seguir dois modelos sem as variáveis mancha e infilt: um com a variável idade e outro sem.
# Com a variável IDADE
modelr3 = glm(trat ~ idade + leuc + tmax, family = binomial)
summary(modelr3)##
## Call:
## glm(formula = trat ~ idade + leuc + tmax, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.76104 -0.68683 -0.09747 0.67388 2.16510
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 87.38804 35.45816 2.465 0.01372 *
## idade -0.05850 0.02558 -2.287 0.02218 *
## leuc 0.38493 0.12152 3.168 0.00154 **
## tmax -0.08897 0.03607 -2.467 0.01363 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 70.524 on 50 degrees of freedom
## Residual deviance: 43.265 on 47 degrees of freedom
## AIC: 51.265
##
## Number of Fisher Scoring iterations: 6
# Sem a variável IDADE
modelr4 = glm(trat ~ leuc + tmax, family = binomial)
summary(modelr4)##
## Call:
## glm(formula = trat ~ leuc + tmax, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2508 -0.9029 -0.1332 0.7621 1.7328
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 70.82281 30.70009 2.307 0.02106 *
## leuc 0.33631 0.10422 3.227 0.00125 **
## tmax -0.07465 0.03144 -2.374 0.01758 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 70.524 on 50 degrees of freedom
## Residual deviance: 49.645 on 48 degrees of freedom
## AIC: 55.645
##
## Number of Fisher Scoring iterations: 5
Em ambos os modelos, a 5% de significâncias, todas as variáveis são significativas.
1 - pchisq(modelr3$deviance, modelr3$df.residual)## [1] 0.6280164
1 - pchisq(modelr4$deviance, modelr4$df.residual)## [1] 0.4075115
Ao nível de 5% de significância, não rejeitamos a hipótese de que os modelos sejam adequados.
fit.model = modelr3
source("envel_bino.txt")fit.model = modelr4
source("envel_bino.txt")Pelos gráficos de envelope, temos que os dois modelos ficaram bem ajustados. Então, pelo princípio da parcimônia, optaremos pelo modelo sem a variável idade.
HosmerLemeshowTest(fit=fitted(modelr4), trat, ngr=10, X=cbind(leuc, tmax))$gof##
## le Cessie-van Houwelingen-Copas-Hosmer global goodness of fit test
##
## data: fitted(modelr4) and trat
## z = 0.006054, p-value = 0.9952
Pelo teste de Homser-Lemeshow, ao nível de 5% de significância, não rejeitamos a hipótese de bom ajuste do modelo.
par(mfrow=c(2,2))
# ALAVANCAGEM
cut = min(sort(hatvalues(modelr4), decreasing = TRUE)[1:3])
plot(hatvalues(modelr4), pch = 16,
xlab = "observação", ylab = "medida h", main = "Alavancagem")
text(hatvalues(modelr4) * as.integer(hatvalues(modelr4) >= cut),
cex=0.8, p=1, offset=0.4)
# DISTÂNCIA DE COOK
plot(modelr4, which=4, lwd=3, main = "Distância de Cook", caption = F)
# RESÍDUO COMPONENTE DO DESVIO POR ÍNDICE
plot(residuals(modelr4), pch = 16,
ylab = "res. compon. do desvio", xlab = "índice")
abline(h = c(-2,2), col="red")Pelo gráfico de resíduo componente do desvio, verificamos que os pontos estão bem distribuídos aleatoriamente no intervalo (-2,2) e sem pontos aberrantes para muito além dele. Nos gráficos de alavancagem e distância de Cook, destacam-se os pontos 2, 32, 40 e 42. Então, a seguir vamos verificar a se há influência desses pontos sobre a inferência do modelo.
summary(influence.measures(modelr4))## Potentially influential observations of
## glm(formula = trat ~ leuc + tmax, family = binomial) :
##
## dfb.1_ dfb.leuc dfb.tmax dffit cov.r cook.d hat
## 2 -0.72 0.02 0.71 0.92_* 1.16 0.25 0.23_*
## 32 0.27 -0.44 -0.26 -0.80_* 1.04 0.22 0.16
## 40 -0.12 0.20 0.12 0.36 1.22_* 0.03 0.16
## 48 -0.35 -0.49 0.36 -0.55 0.80_* 0.22 0.05
Observação 2
modelr4_sem2 = glm(trat ~ leuc + tmax, subset = -2, family = binomial)
impacto(modelr4, modelr4_sem2)[3:5]## impacto pval_com_i pval_sem_i
## (Intercept) -31.779726 0.021058970 0.010676263
## leuc -1.303485 0.001251647 0.001803289
## tmax -30.540927 0.017582847 0.009263493
A 5% de significância, removendo-se a observação 2, não há mudança inferencial no modelo.
Observação 32
modelr4_sem32 = glm(trat ~ leuc + tmax, subset = -32, family = binomial)
impacto(modelr4, modelr4_sem32)[3:5]## impacto pval_com_i pval_sem_i
## (Intercept) 11.05764 0.021058970 0.0449955910
## leuc -15.84809 0.001251647 0.0009397309
## tmax 10.00241 0.017582847 0.0367816803
A 5% de significância, removendo-se a observação 32, não há mudança inferencial no modelo.
Observação 40
modelr4_sem40 = glm(trat ~ leuc + tmax, subset = -40, family = binomial)
impacto(modelr4, modelr4_sem40)[3:5]## impacto pval_com_i pval_sem_i
## (Intercept) -4.298702 0.021058970 0.017876120
## leuc 4.587500 0.001251647 0.002136859
## tmax -3.948505 0.017582847 0.014971725
A 5% de significância, removendo-se a observação 40, não há mudança inferencial no modelo.
Observação 48
modelr4_sem48 = glm(trat ~ leuc + tmax, subset = -48, family = binomial)
impacto(modelr4, modelr4_sem48)[3:5]## impacto pval_com_i pval_sem_i
## (Intercept) -28.40658 0.021058970 0.012766838
## leuc -28.35986 0.001251647 0.001124342
## tmax -28.19501 0.017582847 0.010869361
A 5% de significância, removendo-se a observação 48, não há mudança inferencial no modelo.
Com essa análise, chegamos ao seguinte modelo para explicar a probabilidade de eficiência do tratamento:
\[ \log \left\{ \frac{\pi(\mathbf{x})}{1-\pi(\mathbf{x})} \right\} = 70.82 + 0.34 \cdot leuc - 0.07 \cdot tmax \; . \]
Como ilustração, ao fixarmos a temperatura máxima e considerarmos um aumento de 1% no número de células leucêmicas, a razão de chances de tratamento satisfatório, denotada por \(\psi_{CL}\), é estimada em
\[ \psi_I = e^{0.34} = 1.40 \; . \]
Isto significa que, sob uma mesma temperatura máxima, o aumento de 1% no número de células leucêmicas representa uma aumento de 40% na chance do tratamento ser satisfatório.
Exemplo disponível em: Paula, G. A. (2013). Modelos de Regressão com Apoio Computacional. São Paulo, SP: IME-USP. (Exercício 6, página 345). Dados disponíveis em: https://www.ime.usp.br/~giapaula/textoregressao.htm, arquivo canc1.dat.
Os dados a seguir são provenientes de um estudo de seguimento para estudar a associação entre a taxa anual de câncer nasal em trabalhadores de uma refinaria de níquel do País de Gales e algumas variáveis explicativas: idade no primeiro emprego (4 níveis), ano do primeiro emprego (4 níveis) e tempo decorrido desde o primeiro emprego (5 níveis). São também apresentados o número de casos de câncer nasal e o total de pessoas-anos para cada combinação desses três fatores. Nesse trabalho foi proposto um modelo log-linear com resposta de Poisson, sendo o número de casos de câncer nasal com offset dado por log(pessoas-anos).
As variáveis e suas classes:
idade: 1 (<20), 2 (20-27), 3 (27.5 - 34.9), 4 (35+), em anos;
ano: 1 (<1910), 2 (1910-1914), 3 (1915-1919), 4 (1920-1924), em anos;
tempo: 1 (0-19), 2 (20-29), 3 (30- 39), 4 (40-49), 5 (50+), em anos.
Denotamos por \(Y_{ijk}\) o número de casos de câncer nasal para o i-ésimo nível de idade no primeiro emprego, para o j-ésimo ano do primeiro emprego e para o k-ésimo tempo decorrido desde o primeiro emprego, \(i,j=1,...,4\) e \(k=1,...,5\). Supomos, então, que
\[ Y_{ijk} \sim Poisson(\lambda_{ijk} t_{ijk}), \]
em que \(\lambda_{ijk}\) é a taxa média de casos de câncer nasal por unidade de tempo para a idade i, ano j e tempo k. O modelo saturado, nesse caso, é dado por
\[ \log(\lambda_{ijk}t_{ijk}) = \alpha + \beta_i + \gamma_j + \delta_k, \]
em que \(\beta_i\) é o efeito da i-ésima classe de idade no primeiro emprego, sendo \(\beta_1 = 1(<20)\); \(\gamma_j\) o efeito da j-ésima classe de ano do primeiro emprego, sendo \(\gamma_1 = 1(<1910)\); \(\delta_k\) o efeito da k-ésima classe de tempo decorrido desde o primeiro emprego, sendo \(\delta_1 = 1(0-19)\).
##
## Call:
## glm(formula = ncasos ~ idade + ano + tempo + offset(log(total)),
## family = poisson)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7941 -0.7074 -0.3282 0.2736 2.4937
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -9.27178 1.31892 -7.030 2.07e-12 ***
## idade2 1.67276 0.75214 2.224 0.02615 *
## idade3 2.48173 0.75909 3.269 0.00108 **
## idade4 3.42817 0.78160 4.386 1.15e-05 ***
## ano2 0.61895 0.37130 1.667 0.09552 .
## ano3 0.05414 0.46808 0.116 0.90791
## ano4 -1.12614 0.45294 -2.486 0.01291 *
## tempo2 1.59815 1.04748 1.526 0.12708
## tempo3 1.75121 1.05521 1.660 0.09700 .
## tempo4 2.35486 1.07010 2.201 0.02776 *
## tempo5 2.81758 1.11828 2.520 0.01175 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 135.681 on 71 degrees of freedom
## Residual deviance: 58.166 on 61 degrees of freedom
## AIC: 153.39
##
## Number of Fisher Scoring iterations: 6
Pelos gráficos acima, percebemos que as observações 62 e 63 se destacam nos gráficos de alavancagem e distância de Cook, caracterizando-se como pontos de alavanca e possíveis pontos de influência. A observação 52 se afasta um pouco do limite superior do intervalo de confiança para os resíduos componentes do desvio e será considerada como um ponto aberrante. A seguir, estudaremos a influência desses três pontos sobre o ajuste do modelo.
Observação 52:
## impacto pval_com_i pval_sem_i
## (Intercept) 1.5128475 0.000000000002067477 0.000000000004356923
## idade2 3.6554899 0.026148616816039749 0.032007215605215911
## idade3 5.3753065 0.001078024818957924 0.001991338461317067
## idade4 2.4766547 0.000011541390822584 0.000018039837865180
## ano2 3.4807475 0.095523831910207815 0.108029138724031110
## ano3 272.8270147 0.907913016477262413 0.846449314051831769
## ano4 -1.8553932 0.012908715898225178 0.011467766786010202
## tempo2 0.2290205 0.127083241095111599 0.128300569536079490
## tempo3 0.6315893 0.096997436979855758 0.099518186750867349
## tempo4 0.9039952 0.027764429970676213 0.029392120989161854
## tempo5 7.7644850 0.011749976792630948 0.021946993370582788
Observação 62
## impacto pval_com_i pval_sem_i
## (Intercept) -0.1251849 0.000000000002067477 0.000000000001868126
## idade2 2.1123292 0.026148616816039749 0.029302768207374805
## idade3 1.7230709 0.001078024818957924 0.001300403343998420
## idade4 -6.9946550 0.000011541390822584 0.000003030112064364
## ano2 -14.8348937 0.095523831910207815 0.054868569500244241
## ano3 218.9254636 0.907913016477262413 0.892509547081333787
## ano4 -7.5608046 0.012908715898225178 0.008225427383831935
## tempo2 -8.8173583 0.127083241095111599 0.096874999883543325
## tempo3 2.3170463 0.096997436979855758 0.105239733814455366
## tempo4 -0.2317686 0.027764429970676213 0.027524866857942212
## tempo5 -1.1706799 0.011749976792630948 0.010909475879866461
Observação 63:
## impacto pval_com_i pval_sem_i
## (Intercept) 1.5520620 0.000000000002067477 0.000000000004898349
## idade2 -1.5614680 0.026148616816039749 0.024150242028928428
## idade3 -0.6671066 0.001078024818957924 0.001021541390270312
## idade4 11.1067183 0.000011541390822584 0.000167202020944200
## ano2 28.7232661 0.095523831910207815 0.251882602591362581
## ano3 -67.1709199 0.907913016477262413 0.847118892595951078
## ano4 -0.4606752 0.012908715898225178 0.012672163470177361
## tempo2 -4.4846202 0.127083241095111599 0.110680742091644552
## tempo3 16.0756078 0.096997436979855758 0.168525091205826927
## tempo4 2.8124054 0.027764429970676213 0.032349724301361418
## tempo5 3.5454736 0.011749976792630948 0.014984580367553775
Ao nível de 5% de significância, nenhuma das remoções testadas acima ocasionou mudanças inferenciais nas estimativas dos parâmetros.
Com essa análise, confirmamos o modelo saturado para explicar a taxa anual de casos de câncer nasal.
\[ \log(\lambda_{ijk}t_{ijk}) = \alpha + \beta_i + \gamma_j + \delta_k \]
resumo = summary(model)
table = data.frame(Efeito = c("Constante","20-27","27.5-34.9","35+",
"1910-1914","1915-1919","1920-1924",
"20-29","30-39","40-49","50+"),
Parametros = c("alfa",
"beta2","beta3","beta4",
"gama2","gama3","gama4",
"delta2","delta3","delta4","delta5"),
Estimativa = resumo$coefficients[,1],
E.Padrao = resumo$coefficients[,2])
knitr::kable(table, "pipe")| Efeito | Parametros | Estimativa | E.Padrao | |
|---|---|---|---|---|
| (Intercept) | Constante | alfa | -9.2717833 | 1.3189152 |
| idade2 | 20-27 | beta2 | 1.6727561 | 0.7521394 |
| idade3 | 27.5-34.9 | beta3 | 2.4817313 | 0.7590948 |
| idade4 | 35+ | beta4 | 3.4281689 | 0.7816028 |
| ano2 | 1910-1914 | gama2 | 0.6189461 | 0.3713040 |
| ano3 | 1915-1919 | gama3 | 0.0541431 | 0.4680769 |
| ano4 | 1920-1924 | gama4 | -1.1261375 | 0.4529412 |
| tempo2 | 20-29 | delta2 | 1.5981513 | 1.0474835 |
| tempo3 | 30-39 | delta3 | 1.7512094 | 1.0552076 |
| tempo4 | 40-49 | delta4 | 2.3548632 | 1.0701010 |
| tempo5 | 50+ | delta5 | 2.8175829 | 1.1182812 |
Notamos, então, que as estimativas acima são significativamente diferentes de 0 e que há fortes indícios de que a taxa anual de câncer nasal cresce exponencialmente com o aumento da idade no primeiro emprego e o tempo decorrido desde o primeiro emprego, mas decresce exponencialmente com o aumento de anos do primeiro emprego.